function skeleton = getSkelImage(bwImage, numEndpoints)
 %GETSKELIMAGE - Returns a binary image of the skeleton without spurs
 
 %This algorithm assumes that the longest spurs contain the endpoints of
 %the worm, and that the object has the specified number of endpoints.
 
 skeleton = bwmorph(bwImage, 'skel', inf);
 
 %% Pre-Process the skeleton with spurs
 
 %Remove a pixel from 3x2 blocks that sometimes form when more than one
 %branch spur comes from branchpoints that are adjacent.
 skeleton = filter3by2Blocks(skeleton);
 
 %Remove small loops from teh end of spurs that are sometimes formed when
 %the boundary of the region is very rough and has hangng pixels. 
 skeleton = fixLoops(skeleton,5);
 
 %% Iteratively remove pixels from the end of spurs
 
 %Get a n x 2 matrix of endpoint pixels
 x = bwmorph(skeleton, 'endpoints');
 [row ,col ] = find(x);
 endpoints = [row col];
 
 %Remove pixels form each remaining spur until there are only the specified
 %number of endpoints left
 if size(endpoints,1) >numEndpoints
     
     %Save the pixels belonging to each branch in a cell array with each
     %cell representing a branch
     %Initalize the branches with the first set of pixels
     branches = cell(size(endpoints,1),1 );
     for i = 1:size(endpoints,1);
        branches{i,1} = endpoints(i,:);
     end

     %Remove the pixels from the skeleton that were added to the branch
     index = sub2ind(size(bwImage),endpoints(:,1),endpoints(:,2));  
     skeleton(index) = 0;

     %Get the next set of endpoints
     x = bwmorph(skeleton, 'endpoints');
     [row, col ] = find( x);
     endpoints = [row col ] ;
     
     if size(endpoints,1) >numEndpoints
         
        %Flag to run the loop one extra time to make sure that two branches
        %have more pixels than the rest.
        extrarun = 0;
        
        while size(endpoints,1) > numEndpoints || extrarun == 1

             %Add new endpoints to their repsective branch, but need to
             %make sure that a branch within a branch does not result in
             %two endpoints being added to one branch on a single iteration
             branchMatch = zeros(size(branches,1),1);
             for i = 1:size(endpoints,1);
                 for j =1: size(branches,1)
                      aBranch = branches{j};
                      prevPixel = aBranch(size(aBranch,1),: );
                      prevPixelRow = prevPixel(1,1);
                      prevPixelCol = prevPixel(1,2);
                      if    (endpoints(i,1) == prevPixelRow + 1 || endpoints(i,1) == prevPixelRow - 1  || endpoints(i,1) == prevPixelRow )&&...
                            (endpoints(i,2) == prevPixelCol + 1 || endpoints(i,2) == prevPixelCol - 1 ||endpoints(i,2) == prevPixelCol ) &&...
                            branchMatch(j) == 0;

                          branches{j} = [branches{j}; endpoints(i,:) ];
                          branchMatch(j) = 1;
                          break;
                      end
                 end
             end
             
             %Append zeros to branches that have just ended
             for k =1: size(branches,1)

                 %If no pixels were added to this branch durin the current
                 %iteration
                 if branchMatch(k) == 0 

                     %Get the current branch from the branch list
                     aBranch = branches{k};

                     %Get thelast row of the current branch
                     lastRow  = aBranch(size(aBranch,1),: );

                     %If the last row is not zeros, add a zero
                     if lastRow(1) ~= 0
                        branches{k} = [branches{k}; [0,0] ];
                     end
                 end
             end
                        
             %Delete endpoints from the skeleton image
             index = sub2ind(size(bwImage),endpoints(:,1),endpoints(:,2));
             skeleton(index) = 0;
             
             %Execute the skel morph op once more to get rid of t's
             skeleton = bwmorph(skeleton,'skel',inf);

             %Get the next set of endpoints
             x = bwmorph(skeleton, 'endpoints');
             [row, col ] = find( x);
             endpoints = [row col ];

             %Even if there are only two branches left, run it one extra 
             % time to make sure the head and tail have the most pixels in
             %their branches

             %If it already ran one extra tme, stop looping
             if extrarun == 1
                 extrarun = 2; %stop
             end

            %If it needs to run one more time, do an extra loop
            if size(endpoints,1) <= numEndpoints && extrarun == 0
                 extrarun = 1;
            end
             
        end
        
        %Add the head/tail branches back to the skeleton image 
        for i = 1: size(branches,1) 

            %Get the current branch from the branch list
            aBranch = branches{i};

            %Get thelast row of the current branch
            lastRow  = aBranch(size(aBranch,1),: );

            %If the last row is not zeros, add its pixesl back to
            %the skeleton image
            if ~ isequal(lastRow, [0,0])
               index = sub2ind(size(bwImage),aBranch(:,1),aBranch(:,2));  
               skeleton(index) = 1;
            end
        end
         
    end
   
 end
 
 skeleton = bwmorph(skeleton, 'skel');
 
 %Execute the skel spur opeation to remove sngel pixel spurs that sometime 
 %remain.  
  x = bwmorph(skeleton, 'endpoints');
 [row ,col ] = find(x);
 endpoints1 = [row col];
 if size(endpoints1,1) > numEndpoints
    skeleton2 = bwmorph(skeleton,'spur');
    x = bwmorph(skeleton2, 'endpoints');
    [row ,col ] = find(x);
    endpoints2 = [row col];
    dontDelete = [];
    for j = 1:size(endpoints1,1)
        for k = 1:size(endpoints2,1)
            if isNeighbor(endpoints1(j,:), endpoints2(k,:))
                dontDelete = [dontDelete;endpoints1(j,:)];
            end
        end
    end
    [~, loc] = ismember(dontDelete, endpoints1, 'rows');
    toDelete = [];
    for i = 1: size(endpoints1,1)
        if ~ismember(i,loc)
            toDelete = [toDelete;endpoints1(i,:)];
        end
    end
    index = sub2ind(size(bwImage),toDelete(:,1),toDelete(:,2));
    skeleton(index) = 0;
     
 skeleton = bwmorph(skeleton, 'skel');


end
    
    
    
    
